%% master_counterfactuals_CES_idiosyncraticprefshock.m

clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM
%
% Date: July 2019
%
% Purpose: This file runs the counterfactual exercises after the
% main_baseline.m has been run. Calibrate idiosyncratic productivity shocks
% that are equivalent in size to firms' exposure to the foreign shock
% 
% Based on one file: Step_4_Algorithm_counterfactual.m, which uses new
% baseline input data, MAT/Step_4_setup_`alphaT''rhoT'_`shockT'.mat
% as initial level variables, and then shocks the system and computes hat
% variables for wages, prices, exports, etc.
%
% NEW: for CES setup. Have streamline file that has to be called up and
%                     various variables that have to be used.
% Output: MAT/Step_4_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat,
%         MAT/Step_5_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%         MAT/Elasticity_output_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Loading in the previous data to get some useful variables

alphaT = 4;
rhoT = '3';
lambdaT = '1.000001';
etaT = '1.000001';
lambda = 1.000001;
eta = 1.000001;
sigmaT='1.5';



fname = strcat('MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat');
load(fname,'FI_J', 'firmsecD_sorted', 'start_sorted');
load('MAT/WIOD.mat','pi_c_mnj','countrysecD','france','J','N','start','sum_by_country_dummy','gamma','alpha_WIOD_nj');


rho = 3*ones(J,1);
sigma = 1.5*ones(J,1);

%% SET PARAMETERS VALUES

% Frisch elasticity for GHH preferences

varphi=3;

% share of labour in total final consumption (need to make assumption or collect data)

sL_n = ones(1,N)*0.8; 

% share of profits in total final consumption (should have this data)

sPi_n= ones(1,N)*0.15; 

% share of trade deficit in total final consumption (should have this data)

sD_n= ones(1,N)*0.05; 

%% Baseline: no shocks

% Pi_hat: Change in profits per country
    
Pi_hat_n=ones(1,N);
    
% D_hat: Change in trade deficit
    
D_hat_n = ones(1,N);
    
% Taste shock - for France only
    
Xi_hat_mnjf=ones(FI_J,N);
    
%Taste shock - for sectors only (countries except France)
    
Xi_hat_mnj= ones(N*J,N);
    
% Productivity shock, a is firm specific - France only
    
a_hat_f= ones(FI_J,1);
    
% Productivity shock, a is for sectors - countries except France
    
a_hat = ones(J*N,1);

%% creates an index that captures each sector-decile set of firms
eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',sigmaT,'/Step_5_output_alpha',num2str(alphaT),...
    '_rho',rhoT,'_WORLDprod_p10','_approx.mat')]);

ps=1;
PS=[];
for j=1:J
    perc_sec = zeros(size(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)));
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),99);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 1;
    p = prctile(full(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)),90);
    perc_sec(VA_fj_0(start_sorted(j):start_sorted(j+1)-1,1)>=p & perc_sec ==0) = 9;
    perc_sec(perc_sec ==0) = 90;
    PS = [PS;perc_sec];
end

shockT = 'DOMidioprod';
for sim=[1 9 90]
    
    Pi_hat_n=ones(1,N);
    D_hat_n = ones(1,N);
    Xi_hat_mnjf=ones(FI_J,N);
    Xi_hat_mnj= ones(N*J,N);
    
    a_hat = ones(J*N,1);
    a_hat_f=ones(FI_J,1);
    a_hat_f(PS==sim,1)=.9091;
 
        %% Running the Algorithm

        disp('BEGINNING COUNTERFACTUAL EXERCISE')

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
        % Step_4_Algorithm: same algorithm as Step_4_Algorithm_baseline.m, but now
        % we are able to feed in the model the calculated t+1 varialbes as t,
        % starting at the new equilibrium with no wedges
        eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
    '_rho',rhoT,'_approx.mat')]); 

         Step_4_Algorithm_counterfactual_CES_approx
% 
%         % Save the counterfactual output so it can be used for post analysis
% 
         outputtable = [' X_mnj X_hat_mnj pi_mnjf1 pi_mnj1 ' ...
                  'pi_c_mnj1 pi_c_nj1 pi_M1 pi_l1 pi_M_f1 pi_l_f1 '...
                  'PC_hat_n1 PC_n P_hat_n P_hat_nj P_hat_mnj w_hat ' ...
                  'a_hat a_hat_f Xi_hat_mnj Xi_hat_mnjf '...
                  'final_goods1 intermediates1'] ;
%         %  
         eval([strcat('save MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
                  '_rho',rhoT,'_',shockT,'_approx_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]);   
%         %     
         display('STEP 4 COUNTERFACTUAL COMPLETED SUCCESFULLY')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


            % Step_5_counterfactual_output: Allows the user to easily analyse the
            %      counterfactual exercises produced before. Runs the CounterOutput.m
            %      function after specifying alpha type (labor share), gamma
            %      calibration and the shock scheme to be analysed.

            %% Load data

            
            % Time 0 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
     '_rho',rhoT,'_approx.mat')]);

            % Time 1 data

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
     '_rho',rhoT,'_',shockT,'_approx_Percentile',num2str(sim),'.mat')]);
 
        Step_5_counterfactual_output_CES_approx

            % Save the counterfactual output so it can be used for post analysis
        outputtable = [' VA_hat_n VA_hat_mj VA_hat_mnj VA_hat_fj VA_hat_fnj '...
          'RGDP_def_hat_n RGDP_def_hat_fj RGDP_cpi_hat_n RGDP_cpi_hat_fj '...
          'RGDP_doubledef_hat_n RGDP_doubledef_hat_fj '...
          'real_wage_cpi_hat real_wage_def_hat real_wage_doubledef_hat '...
          'real_wage_income_cpi_hat_n real_wage_income_def_hat_n real_wage_income_doubledef_hat_n '...
          'TFP_hat_def_n TFP_hat_cpi_n TFP_hat_doubledef_n TFP_hat_cpi_france_j '...
          'imports_over_GDP_all imports_over_GDP_france_sector '...
          'exports_over_GDP_all exports_over_GDP_france_sector '...
          'GDP_deflator_n DoubleDeflation_deflator_n P_hat_n VA_fj_0 VA_mj_0'] ;
        %  
        sigmaT = '1.5';
             eval([strcat('save MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]); 

        %     
             display('STEP 5 COMPLETED SUCCESFULLY')

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % Step_6_Elastisticity_output: perform covariance decomposition at firm and
            %        sector level

            % Load files: only needed if run separately from Step 5

            eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Percentile',num2str(sim),'.mat;')]); 

            % Shock size

           delta = 0.1;

           Step_6_Elasticity_output_CES

       outputtable = [' e_Y_FRf_cpi m_e_f_cpi cov_e_f_cpi check_f_cpi '...
           'e_Y_FRs_cpi m_e_mjFRs_cpi cov_e_mjFRs_cpi check_mjFRs_cpi '...
           'e_Y_s_cpi m_e_mj_cpi cov_e_mj_cpi check_mj_cpi '...
           'e_Y_FRf_def m_e_f_def cov_e_f_def check_f_def '...
           'e_Y_FRs_def m_e_mjFRs_def cov_e_mjFRs_def check_mjFRs_def '...
           'e_Y_s_def m_e_mj_def cov_e_mj_def check_mj_def '...
           'e_Y_FRf_doubledef m_e_f_doubledef cov_e_f_doubledef check_f_doubledef '...
           'e_Y_FRs_doubledef m_e_mjFRs_doubledef cov_e_mjFRs_doubledef check_mjFRs_doubledef '...
           'e_Y_s_doubledef m_e_mj_doubledef cov_e_mj_doubledef check_mj_doubledef'] ;

       eval([strcat('save MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
                 '_rho',rhoT,'_',shockT,'_approx_Percentile',num2str(sim),'.mat',outputtable,' -v7.3;')]); 

             display('STEP 6 COMPLETED SUCCESFULLY')

        disp('COUNTERFACTUAL EXERCISE COMPLETED SUCCESFULLY')

end

 